Feldman-Cousins Confidence Levels — Toy MC Method 
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In particle physics, the likelihood ratio ordering principle is frequently used to determine con- 
fidence regions. This method has statistical properties that are superior to that of other confi- 
dence regions. But it often requires intensive computations involving thousands of toy Monte Carlo 
datasets. The original paper by Feldman and Cousins contains a recipe to perform the toy MC 
computation. In this note, we explain their recipe in a more algorithmic way, show its connection 
to 1 — CL plots, and apply it to simple Gaussian situations with boundaries. 
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I. LIKELIHOOD RATIO ORDERING 

Feldman and Cousins have developed a method 1] to 
compute confidence intervals that has well-defined fre- 
quentist coverage. Their method has certain advantages 
over the classical approach to compute confidence levels: 
it avoids the so-called "flip-fiopping" because it offers a 
natural transition between upper (or lower) limits and 
two-sided confidence intervals. It usually gives a non- 
empty interval even in cases, in which an observable is 
measured to be far in a forbidden region-a situation in 
which classical constructions may lead to empty intervals, 
that don't contain the true value with absolute certainty. 
And the method effectively decouples the goodness-of-fit 
confidence for that for an estimated parameter (cf. the 
instructive example in Ref. Q). 

In general, frequentist confidence intervals [/ii,/i2] of 
confidence level (CL) a are computed by constructing 
the confidence belt (Neyman construction, cf. Ref. [3]). 
This belt consists of the conjunction of the intervals 
[xi{fi), X2{fJ.)], that is for each value of /i one finds an 
interval by integrating the probability density function 
P(a;|/i) such that 



(1) 



The belt has the property, that as long as each of the 
intervals [a::i,a;2] fulfill Eq. [U each vertical (that is in fi 
direction) intersection gives an interval of coverage a. 
Thus, after the belt is constructed, one determines the 
confidence interval [fii , /i2] by intersecting the belt at the 
measured value of x. 

But Eq. [T] doesn't determine xi and X2 uniquely. In 
order to solve uniquely for xi and X2 one also needs to 
specify an ordering principle that determines in which or- 
der the dx shall be included into the interval. The classi- 
cal intervals are obtained by including the dx in order of 
their probability. This yields the additional requirement 
that 



P(xi|/i)=P(x2|M) 



(2) 



Feldman and Cousins define an alternative ordering 
principle. They include the dx in oder of decreasing like- 



lihood ratio, 



Rix,fj.) 



Pix\fihcst) 



(3) 



Here, /Xbcst is the value of /i that maximizes the likelihood 
P{x\fj,) and at the same time is inside the allowed region. 
Note the distinction between probability and likelihood. 
Although numerically identical, one refers to P{x\fi) as 
probability density only if interpreted as a function of x. 
If interpreted as a function of the true parameter fi, it 
is called likelihood, because in the frequentist interpreta- 
tion the true parameter is what it is and it doesn't make 
sense to assign a probability to it. In practice, using the 
likelihood ratio ordering means that for a given value of 
fi, one finds the interval [xi, X2] such that 



Rixi) = Rix2) 



(4) 



and Eq. [T] hold. It is worth noting that in the case of 
discrete probability densities such as the Poisson distri- 
bution, the equalities of Eqns.[T]|3]in general can't be met 
exactly. 



II. TOY MC METHOD BY F.-C. 

In their original publication [ij Feldman and Cousins 
offer a method to solve Eqns. [T] and S] numerically. The 
method is based on toy Monte Carlo (MC) experiments, 
which are datasets drawn from the assumed p.d.f. They 
apply the method to a non-trivial example from neutrino 
physics. Here we shall use it to compute the confidence 
belt for the trivial example of a unit Gaussian bound to 
positive values, /i > 0. 

In practice it is convenient to take the logarithm of the 
likelihood ratio 



-21ni?= Ax' = x'-xL 



(5) 
(6) 



because one often prefers to minimize — 2 In P rather than 
to maximize P. Due to the monotonic nature of the 
logarithm, neither the position of the extrema nor the 
ordering principle change. 
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The algorithm proposed by F.-C. computes the inter- 
val [xi, X2] for each value of the true parameter (the mean 
/i of the Gaussian) in a sufficiently fine grid. The confi- 
dence belt is the conjunction of all such intervals. The 
algorithm is described in turn: 

1. At the considered value of the true parameter /ip, 
generate a toy experiment by drawing a value Xtoy 
from the p.d.f. That is, draw from the unit Gaus- 
sian G{a — 1, /i — /io)- 

2. Compute ^Xtoy the toy experiment, follow- 
ing Eq. [6] For x^(x,/i), x is set to a;toy and 
fi is set to ^0- For x'^ , fJ-bcst) , x = Xtoy and 
/^bcst — max(0,xtoy) implementing the boundary 
at zero. 

3. Find the value Axl, such that a of the toy experi- 
ments have Axtoy < ^xl- 

4. The interval [a;i, X2] is given by all values of x such 
that Ax^(a;,Aio) < Axl 

Note that the general case quickly gets more expensive 
than the sketched example. If there are nuisance param- 
eters present, the likelihood needs to be minimized in 
their respect when calculating and Xbost- Also, the 
fact that in step 2 we can set /ibost — a^toy if away from 
the boundary is a consequence of the x -ir^ fi symmetry 
of a Gaussian. In general, fj-hest needs to be obtained 
from a minimization, that respects the boundaries of the 
problem. 

The algorithm is illustrated in Fig. [T] It displays all 
toy experiments generated in step 1, and the curve built 
by their Axtoy values. At values for the measured x 
greater than zero, it follows a parabola, Ax^ = {x — . 
At values below zero the effect of the boundary becomes 
visible when the curve becomes linear, Ax^ — (x — /i)^ — 
x^ = — 2^x. In step 3, an interval is selected such 
that a fraction a of the toy experiments fall inside. This 
satisfies Eq. [TJ The value Axc for which this is the case 
is depicted in Fig. [1] as a horizontal line. Finally in step 4 
the interval boundaries xi and X2 are determined as the 
intersection of the Axtoy curve and the Axc line. This 
satisfies Eq. 01 The full resulting confidence belt is shown 
in Fig. H 

For illustration, we compute confidence belts for three 
additional similar situations: 90% CL without bound- 
aries, 90% CL with a two-sided boundary < /i < 3, 
and 68% CL with a boundary at > 0. The previous 
example of Fig. [5] is repeated for comparison. 

III. 1-CL PLOTS 

Confidence belts might not always be the most conve- 
nient way to display information about the statistical pre- 
cision of a measurement. For instance, if we performed 
not only one measurement of the above Gaussian mean. 
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Measured Mean x 

FIG. 1: Example set of 5000 toy experiments drawn from 
G(/i — 0.5, a — 1) (histogram), together with the likelihood 
ratio following Eq. [6] computed in presence of a boundary 
fi > (solid line) , and the value Axc giving a = 90% (dashed 
line). 
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FIG. 2: Confidence belt for a unit Gaussian with a mean 
bound to be positive, a = 90%. 

the belt becomes a multi-dimensional object. Instead, 
one often choses to show "1 — CL" plots for a parameter 
of interest. An example is given in Fig. SI These plots 
display confidence regions for all possible confidence lev- 
els, evaluated for one specific (the observed) measured 
value xq. To read off the, say, 68% CL interval [/ii,/i2] 
one draws a horizontal line at 1 — CL = 1 — 0.68 — 0.32. 
The interval is given by the intersection of this line and 
the curve. At the most probable value for the true pa- 
rameter the curves raises to one. 

The 1 — CL plots are constructed using a slight vari- 
ation of the original toy MC method described above. 
Steps 1 and 2 remain unchanged: 

1. At the considered value of the true parameter /ig, 
generate a toy experiment. 

2. Compute Axtoy = Ax^ixtoy, fJ-o) for the toy exper- 
iment. 

3. Compute Ax^^ta = ^X^i^o, ^J'o) for the measured 
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FIG. 3: Additional example confidence belts for a unit Gaus- 
sian. 



leading to 

= Ax'(xo,A*)<Ax^} . (10) 

Due to the Neyman construction of the confidence belt 
both horizontal and vertical intersections of the belt give 
intervals of coverage ao . Thus it is also 

h = {x\AxHx,f^o)<^xl} , (11) 

which is just the definition of step 4 of the original algo- 
rithm. 

Two example physics analyses where this method is 
used can be found in Refs. [1, 

IV. TWO-DIMENSIONAL CL CONTOURS 




FIG. 4: Plot of 1 — CL of a unit Gaussian for a measured 
mean of a; = 1.4, both without boundaries (dashed) and for 
a boundary at /i > (solid). The horizontal lines indicate 
a = 68% and a = 90%. 



value. 

4. The 1 — CL (= 1 — a(^o)) value is given by the 
fraction of toy experiments that have a larger Ax^ 
than the measured value: 



1 - a(Mo) 



^^(AxLta 



< Ax?oy) 



toy 



(7) 



The full 1 — CL curve is then obtained by scanning for 
each possible value jio- The equivalence to the original 
algorithm can be seen as follows. The confidence interval 
of coverage ao for the true parameter fj, computed by the 
1 — CL algorithm is given by 



{/i I 1 - a{fio) > 1 - ao} 



(8) 



In step 3 of the original algorithm we have defined A^c 
such that ao of the toy experiments have smaller Ax^ 
values, thus we can replace both sides of the inequality, 



N{Ax^{xo,fi) < Axly) 



> 



N{Axl < AxlJ 



toy 



(9) 



It is straight-forward to apply the toy MC method to 
a two-dimensional situation. If there are two, possibly 
correlated true parameters that shall be estimated from 
the data, the 1 — CL plots become three-dimensional. 
From those, CL contours are readily obtained by slicing 
at the desired confidence level 1 — ao. 

We use the toy MC method to obtain CL contours 
for the means /ii and fi2 of a two-dimensional Gaussian. 
We chose the corresponding widths ai^2 different from 
each other, and also a non-zero correlation. The resulting 
contours are tilted ellipses. They are shown in Fig. [5] 
(solid lines). The innermost ellipsis has the property, 
that its projections onto the /ii^2 axes are exactly ±(Ti^2- 
Its confidence level is ao = 39.3%. 

We also consider the effect of boundaries for the true 
parameters /ii,2- In Fig. [5] (dashed lines), we chose the 
allowed region to coincide with the displayed region. The 
resulting contours all cover a smaller area compared to 
the unconstrained case: The boundaries seem to "de- 
flect" the contour, resulting in "squeezed ellipses" . 



V. PROB METHOD 

In many cases, the likelihood function assumes a Gaus- 
sian shape. This is likely the case, when the parameter 
of interest is itself an estimator that is based on a suffi- 
ciently large data sample. Then, Ax^ follows a x^ dis- 
tribution with one degree of freedom by definition-which 
also justifies the use of the symbol in the above. It is 
therefore possible to compute 1 — CL from the cumula- 
tive distribution function of the x^ distribution (compare 
to Eq.ED: 



1 - CL 



(t) dt 



1 



\/2r(i/2) 



toy 



PROB(AxLta,l) 



(12) 

(13) 
(14) 
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FIG. 5: Confidence regions for the mean of a two-dimensional 
Gaussian with ui = 0.4, 02 = 0.6, p = 70%, for a measured 
mean of {xi,X2) ~ (—0.2, 0.2) (data point) with coverage a — 
39.3%, 68.3%, 90%, 95.4% (solid lines from inner to outer). 
The dashed lines show the corresponding regions in presence 
of the boundaries \fj,i\ < 1 and |/i2| < 1. 



With this approach there is no need for toy MC. The al- 
gorithm in Section lnll iust reduces to steps 3 and 4, where 
in the latter the toy MC ratio is replaced by Eq. [TH The 
Prob function is available within the Root framework 
as TMath: :Prob(). 

It should be stated clearly that this method only ap- 
plies in the Gaussian regime. For instance, it does not 
handle the presence of boundaries in a satisfactory way. 
Applied to the situation of Fig. |4] it produces a curve 



that coincides with the unconstrained F.-C. curve every- 
where: When reaching the boundary fi — 0, it suddenly 
drops off to zero. In this situation, correct coverage is 
not guaranteed. 

The clear advantage of the Prob method is that the 
required amount of computation is very reasonable, as no 
toy MC generation is involved. And often it gives a very 
good idea of the F.-C. confidence contours, so it may be 
instructive to try it before a F.-C. analysis is performed. 



VI. CONCLUSION 

In their original publication, Feldman and Cousins in- 
clude a recipe to perform a toy MC computation of their 
confidence regions. In this note, we explain it in an al- 
gorithmic way, that is suited to direct implementation 
into existing fitter frameworks. We show the connection 
of the toy MC method to the frequentist confidence belt 
and to 1 — CL plots. We illustrate the method using sim- 
ple one- and two-dimensional Gaussian situations with 
and without boundaries on their true parameters. 
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